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Abstract 

We present an efficient second order accurate scheme to treat stiff source terms 
within the framework of higher order Godunov's methods. We empioy Duhamel's 
formula to devise a modified predictor step which accounts for the effects of stiff 
source terms on the conservative fluxes and recovers the correct isothermal behavior 
in the limit of an infinite cooling/reaction rate. Source term effects on the conserva- 
tive quantities are fully accounted for by means of a one-step, second order accurate 
semi-implicit corrector scheme based on the deferred correction method of Dutt et. 
al. We demonstrate the accurate, stable and convergent results of the proposed 
method through a set of benchmark problems for a variety of stiffness conditions 
and source types. 
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1 Introduction 



We wish to solve the following system of partial differential equations describ- 
ing a hydrodynamic flow with a stiff (energy) source term 

where D is the dimensionality of the problem, U, F(U), S(U) are the con- 
servative variables, the conservative fluxes and the source term respectively, 
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In the above equations, p is the density, u d the velocity in the d direction, 
E = e + Y, d =i u d /2 is the total specific energy with, e, the specific internal 
energy. A(e,p) is the term describing the source of specific internal energy. 

In the following we consider the case of a stiff source term corresponding to 
an endothermic process, such as occurs in radiative losses. In addition, we 
restrict our analysis to source types that, at least near equilibrium, behave as 
a relaxation law. In the stiff case, the characteristic relaxation time scale for 
S may be much smaller than the CFL time step for the hydrodynamic waves. 
For that reason, we would like to use a semi-implicit method, treating the stiff 
source term implicitly, while using an explicit method for the hyperbolic terms. 
However, the classical analysis of such fast endothermic processes shows that, 
in the limit as the relaxation time goes to zero, the gas can be described by the 
compressible flow equations with an isothermal equation of state [18] . Pember 
[14] showed that the use of formally second-order accurate semi-implicit meth- 
ods such as Strang splitting, or a second-order Godunov predictor-corrector 
method could lead to a substantial loss of accuracy, due to inconsistencies 
between the the characteristic tracing step without sources and the effective 
limiting isothermal behavior. Such inconsistencies between the flux calcula- 
tion with the limiting isothermal equation of state can lead to dramatic errors 
particularly at sonic points. Pember proposed various approaches to the prob- 
lem based on classical relaxation theory. Roe and Hittinger [15] also addressed 
the issues raised here in relation to Godunov's method with stiff relaxation. 
In their approach they split the equations based on a splitting of state space 
into stiff and non-stiff subspaces of the linearized source term to obtain in the 
stiff limit formulations similar to ours. However, neither Pember nor Roe and 
Hittinger did present a complete method that is second-order accurate in both 
the stiff and non-stiff limits, nor did they discuss the extension to more than 
one dimension. 

The problem of hyperbolic system with stiff relaxation has also been consid- 
ered by other authors in the past mostly for one-dimensional systems and 
within the framework of Runge-Kutta based methods of lines. In particular 
Jin [7] designed second order Runge-Kutta type splitting methods with the 
correct asymptotic limit. Jin and Levermore's [8] developed a semidiscrete 
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high resolution method which, in order to ensure the correct asymptotic be- 
havior, employes a linear combination of the conservative fluxes for the homo- 
geneous (i.e. without the relaxation term) and equilibrium system. The fluxes 
are computed with a higher order Godunov's method and the scheme allows 
for a rapid transition between the stiff and nonstiff regimes. As the authors 
point out, however, the upwind property of the scheme is not strictly guaran- 
teed for all stiffness conditions. Finally, Caflisch, Jin and Russo [2] developed 
a scheme for hyperbolic systems with relaxation that is uniformly accurate for 
various ranges of stiffness conditions (see also Ref. [9,13]). 

The aim of this paper is to build a higher order Godunov's method that pre- 
serves the properties of robustness and accuracy across a variety of stiffness 
conditions thus avoiding the problems described in [14]. In particular, in or- 
der to preserve higher order accuracy, we aim for a semi-implicit method that 
corresponds to a standard second-order Godunov method of the appropriate 
hyperbolic problem for the stiff or non-stiff limits. To this end, we use second- 
order accurate deferred corrections method of a type presented in [5] , to obtain 
a semi-implicit corrector that is a special case of the algorithms described in 
[12], although any implicit L-stable second-order one-step method would be 
acceptable. The main new idea in our work is contained in our treatment of 
the predictor step for computing the hyperbolic fluxes, based on the deriva- 
tion of a local effective dynamics using Duhamel's formula. This leads to an 
explicit predictor step that corresponds to that for a conventional second- 
order Godunov method for Eq. (1) in the limit where the relaxation time is 
comparable to or greater than the hydrodynamic CFL time step; and to a 
second-order Godunov method for the isothermal equations in the limit where 
the relaxation time is much smaller than the hydrodynamic time step. Our 
approach is similar to that used in [17] for obtaining a well-behaved numerical 
method for incompressible viscoelastic flows in both the viscous and elastic 
limits; however, the details there are quite different than those for the present 
setting. 

The paper is organized as follows. In section 2 we describe a second order ac- 
curate, semi-implicit corrector method based on the deferred corrections ideas 
presented in [5,12] to be used for the final source term update. In section 3, 
based on Duhamel's formula, we work out a modified formulation of Godunov's 
predictor step and flux calculation suitable for the case of stiff source terms. In 
section 4 we discuss stability issues for our approach, and Section 5 contains 
the extension of the method to the case in which the source term depends 
both on the gas density as well as the internal energy. In section 6 we test the 
performance of the code and demonstrate the accuracy of the method in var- 
ious stiffness conditions. The paper concludes with section 7 where the main 
results of the paper are summarized. 
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2 Semi-Implicit Predictor-Corrector 



Our time-discretization for the source terms is a single-step, second-order ac- 
curate scheme based on the deferred correction ideas in Dutt, Greengard & 
Rokhlin [5]. Given the system of equations (1) 



dU 
~dt 



— V ■ F + S(U) 



(3) 



we aim for a scheme in which an explicit approach is retained for the non-stiff 
conservative hydrodynamic term, V • F, and a implicit method is employed 
for the stiff part of the equation, S. The particular approach is a special case 
of a more general class of semi-implicit methods by Minion [12] . Consider the 
first order system of ordinary differential equations (ODEs) 



(4) 



Y(t = 0) = Y (5) 

with Y E R n , C : R x R n -> R n . In [5], Eq. (4) is reformulated in terms of its 
equivalent Picard integral equation, to which a deferred corrections algorithm 
is iteratively applied. First an error is estimated according to 



Y + C 
Jo 



r, Y{t) dr - Y(t) , < t < At 



(6) 



where, Y(t), is an initial guess to the solution to be corrected iteratively. 
Then a correction is computed by solving the error equation for the correction 
S(t) ee Y(t) - Y(t) : 



5{t) = J* [C [t, Y(t) + 6(t)] - C [r, Y(t)] } dr + i(t) 
Y(t) = Y(t) + 8(t) , < t < At. 



(7) 



To complete the specification of the method, we need to choose a quadrature 
scheme to replace the integrals in time by sums over a finite number of points. 
The choice of quadrature method in the error calculation (6) and of the number 
of iterations determines the accuracy of the method. However, as noted in 
[5], the rate of convergence of the method is independent of the accuracy of 
the quadrature rule used in the correction calculation (7). In particular, for 
stiff systems, one uses a quadrature rule corresponding to backward Euler, 
replacing the integrand by its linear approximation. In the present case, we 
are only interested in second-order accuracy, so we can use the trapezoidal 
rule for the quadrature rule in the error calculation, and iterate only once. 
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Our semi-implicit method will correspond to solving a collection of ODEs, one 
at each grid point 



(V-F) 



n+i 



(8) 



where we view the time-centered flux divergence as a constant source, whose 
computation using a modified Godunov method is decried below. Following 
[12], we solve the resulting collection of ODEs using the method described 
above. For our initial guess, we use 



U = U Q + (I - AtVuS] 



U 0J 



S(U ) - (V • 



At 



(9) 



where U = U(t ). In the above expression we have used backward Euler to 
estimate the effects of the source term and we have then Taylor expanded 
the implicit part of it into a linear form. This yields a second order accurate 
estimate in the sense that: U — U(t + At) = 0(At 2 ). Based on Eq. (6) the 
error is then estimated as 



e(At) =U + y [S(U) + S(U )] - At (V • F) n+ ^ - U (10) 



where we have used the trapezoidal rule to estimate the integral of the source 
term. The sought correction is obtained in implicit form by applying backward 
Euler to the integral in the correction equation (7) 



5(At) = (I - AtVuS\ )- L e(At) 
U(t + At) = U + 5(At) 



(11) 
(12) 



From Eq. (10)-(11) it is clear that the final solution will have a truncation 
error 0(At 2 ) and global second order accuracy in time. 

Clearly in the non-stiff limit, as the contribution from the term 'AtV^S'l^' 
becomes negligible compared to those from T, the above scheme reduces to 
the usual second order accurate explicit formulation 



U(t + At) = Uq - At (V • F)" + 5 + ^ \S(U) + S(U { 



(13) 



3 Effective Dynamics and a Modified Godunov's Method 



In order to compute the flux divergence (V • F)™ + 2, we use the quasilinear 
form of the equations in primitive variables to extrapolate from cell centers to 
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cell faces 



SW =VuWS(U). 

In order to develop our formulation we will start using W = (p, u, e) T , but will 
switch to the usual set of primitive variables later in Sec. 3.1.1. Hereafter, we 
will denote = S, dropping the superscript. We can also give the evolution 
along the Lagrangian trajectories 

£ + £<-™ 

DW dW 
A L d =A d -u d I,— = — + (u-V)W 



We will derive from the quasilinear form of the equations a new system that 
includes, at least locally in time and state space, the effects of the stiff source 
terms on the hyperbolic structure, and use that quasilinear system to extrap- 
olate from cell centers to faces in a Godunov method. 

We first illustrate the approach for the case of a system of ODE. Consider the 
system of differential equations 

dY 

— = BY + C(t), Y(t ) = Y (14) 
F:R^M n , BeR nxn , C:R^R n . (15) 

The evolution of the rate of change of Y(t), namely 5Y = Y(t) — Yq, is then 
described by d5Y/dt = B5Y + BY + C(t) with SY(0) = 0. According to 
Duhamel's formula, 

5Y{t) = f e^ B \BY Q + C(r)] dr. (16) 
Jo 

When the properties of B lead to a stiff numerical problem, the exponential 
term in the above integral is the one that changes most rapidly, motivating 
the approximation 

5Y(t) w l B (v)[BY + C(0)]f (17) 

where 

X B (7 7 ) = 7 7 - 1 [ V e TB dr. (18) 
Jo 

*^=I B (ri)[BY + C(0)]- (19) 



6 



In what follows, we will take r] = O(At). There are two distinguished limits 
to the effective equation. First, if \\Br)\ \ <C 1, then 

X B (ri) = 7 + 0(77) (20) 

The second is when there is a single eigenmode of B that is stiff relative to 
the time scale defined by r\. Specifically, we assume that for some v G M n 

B — B — \vv T Bv = 0, v T B = (21) 

with 

A?7»l, ||^||«1. (22) 
Here A -1 is the fast time scale that is stiff relative to rj. So we can write 

l B (r ] ) = (I-vv T ) + o( V ,^j (23) 

and 

X B (ri)B = B + oU\). (24) 



A, 

In this case, T B {rj) projects out the stiff dynamics, leaving only processes that 
are resolved on the 0(r]) timescale. 

We can use the effective equation (19) with rj = At to compute a first-order 
accurate predictor step in a second-order accurate predictor-corrector. 

Y cS (At) = Y(0) + l(At)(BY(0) + C(0))At (25) 

Then 

Y cS (At) - Y(At) = 0(At 2 ) if (23) holds 

= (A* a + ^) if (24) holds (26) 

We apply this idea to the dynamics along Lagrangian trajectories. We define 
5W = W[x(t),t] - W[x(t ),t ] = W-W (27) 

and 

D5W 

__ + G = S + S 5W (28) 
d=i UXd 

We have linearized the source term around the value of the state at the begin- 
ning of the Lagrangian trajectory, with S = Vw • S. By applying Duhamel's 
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formula to Eq. (28) we obtain 



SW(t)= / e^- T ^°(-G + S )dr. 

Jto 

Following similar reasoning to the ODE case, we obtain 

DW cS / D 



Dt 



(30) 



/ D 8W\ 
+ \£W^ W J=WS°- (31) 



3.1 Characteristic Analysis 



We will use the quasilinear system (31) with (77) = At/2 to compute the Go- 
dunov predictor step. In the non-stiff limit, this leads to an 0(At 2 ) error in 
the predicted values at cell faces which is sufficient for second-order accuracy 
in the overall method. In order to do that, we need to analyze the hyperbolic 
structure of those equations. Without loss of generality in the following sub- 
sections we still consider the 1-dimensional case. Also, in this section we will 
focus on the case A p = dA/dp = 0, A e = dA/de ^ 0; we will discuss the more 
general case in Sec. 5. With this choice of S , from Eq. (18) we obtain 



( \ 0^ 



^o( A */2) 



1 
0a 



(32) 



where 



1 



A e At 



a = 



< a < 1. 



(33) 



Thus, the presence of a stiff source term leads us to the transformations: 



A = A L + u\ -> A 



eflf 











P \opJ e P \ oe 



V 







a'- 







ul. 



(34) 



/ 



3.1.1 Modified eigenvalues 

Characteristic analysis of the matrix A eS leads to the characteristic equation 



det(A cfr - AI) = (A-u) 



p ( dp" 



dp" 



de 







(35) 
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which admits the the familiar solutions 
Aq = u, A± — u ± 




(36) 



It appears from the above equation that the presence of the source term alters 
the sound speed according to 




P ( dp\ ( dp^ 
a pAde) \dp. 



For a 7-law equation of state, we have 

P = (7 - l)pe 



(37) 



(38) 



Cos = \[a (7 - 1) + 1] 



P 



(39) 



Thus, in the limit of a negligible source term, a — > 1, c e ff — > ('jp/p) 1 ^ 2 and 
the polytropic behavior is recovered. However, in the limit of a stiff source 
term, a — > 0, c e s — > (p/p) 1 ^ 2 and the isothermal regime is approached. This is 
also apparent from the expression for the rate of change of the internal energy 
along Lagrangian trajectories 



De 
Dt 



—a 



p du 
p dx ' 



(40) 



suggesting the limit, de — > as a — > 0. Notice that in our approach we retain 
the polytropic form of the equation of state p = (7 — l)pe, 7 7^ 1, but we 
avoid differentiating it when the presence of source terms must be taken into 
account. Based on Eq. (40) the pressure change is found to be 



Dp 
~D~t 



Dp 

CcS Dt 



~PC 2 cS 



du 
dx 



(41) 



Finally, we note that in general, in D— dimensions, the above analysis ap- 
plies unaltered to the linear operator, Af, for each direction, d, after properly 
transforming u — > Ud, x — > x d - In addition, D — 1 equations are added de- 
scribing the passive transport of momentum components perpendicular to the 
d direction, and the eigenvalue A acquires multiplicity D. 



3.1.2 Modified Eigenvectors 

Given Eq. (41) we can now replace internal energy with pressure and find out 
the expression for the eigenvectors for the usual set of primitive variables. This 
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reads 



(42) 



where in addition to density, velocity and pressure, we have also included the 
specific entropy, s = pp^ 1 (useful, e.g., for the case of hypersonic flows [11]). 
The change in specific entropy is given by 



Ds = /Dp_ c2 D^ 
Dt ' \Dt Dt. 



)=-^(4-^)|-^| (43) 



The linear operator is 



A 



cff 



Op 

o p- 1 

pc 2 cS 

^o <w 1-7 o oy 



ul. 



(44) 



The extra variable 's' results in an additional eigenvalue, X — u, for the oper- 
ator A eS . The set of left and right eigenvectors are given respectively by 



h= fo,--^,-V,0 N 

V 2c cff 2c 2 eS , 



l 2 = 1,0,- — ,0 

V C eff / 

h= fo,^-,-i-,o 



2c eff 2ci s 



(45) 
(46) 
(47) 
(48) 



£sff 
P 

r 2 



\5 c 2p "<) 



T2 



fl\ 




(o) 












; r 3 = 



















r 4 



Ccff 
P 

c cff 



(49) 



\6 c 2p ij 
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3.2 Godunov Predictor in One Dimension 



With the operator A and the sets of left and right eigenvectors that we have 
worked out in the previous section, the Godunov predictor step is carried out 
as usual as follows. 

First the local slopes are defined. In particular at each point left and right 
one-sided slopes as well as cell centered slopes are evaluated and then a final 
choice on the local slope AWi is defined by using van Leer limiter. The upwind, 
time averaged left (— ) and right (+) states at cell interfaces due to fluxes in 
the normal direction, d, are then reconstructed as: 



W h± = W?+ l -(l- ^-Af ) P ± (AW t ) (50) 

where 

P±(W)= ]T (l k -W)-r k . (51) 

±A fe >0 

The source term component is likewise accounted for as 

W l:±4 = W h±4 + ^l So (At/2)S . (52) 



The fluxes at the cell faces F i+ i are computed by solving the Riemann prob- 
lem with left and right states given by (W i+ ,W i+ i_) to obtain W. i 2 and 
computing F i+ i = F (w^i 2 ^- 

To modify this procedure to account for the effective dynamics, we use the 
characteristic analysis of the effective dynamics to perform each of the three 
steps. The projection operator and any limiting in characteristic variables is 
done using the eigenvectors and eigenvalues for the effective dynamics derived 
in Sec. 3.1. Typical approximate Riemann solvers use weak- wave approxima- 
tions to compute the jumps, which only require the linearized jump relations 
provided by the characteristic analysis for the effective dynamics. For the 
case of a polytropic gas, one can use more nonlinear approximate Riemann 
solvers, e.g. two shock approximations, to compute the jump relations, treating 
1 + 01(7 — 1) as an effective polytropic 7. This is done for the results presented 
here. Finally, any entropy fixes required to eliminate rarefaction shocks require 
only the sound speed, for which we again use c e ff. 
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3.3 Extension to More than One Dimension 



For directionally unsplit schemes in D dimensions an additional step is re- 
quired in order to correct the time-averaged left/right states at cell interfaces, 
Wi^±^d in Eq. (52), for the effects of D — 1 fluxes perpendicular to the cell 
interface normal direction. Based on Eq. (31) the effect of the stiff source 
term would be accounted for by carrying out for each additional direction, d, 
a transformation 

A d - l So (At/2)A L4 + u d l = Af (53) 

analogous to that described in Eq. (34). In the method proposed by [4,16] the 
corrections due to transverse fluxes are computed according to a conservative 
scheme. For example in two dimensions 

W hJ , ± , x = W hJ , ± , x - ^VuW (Fl H - F%_i) (54) 

where the input Wij,± tX is computed using a one-dimensional Godunov calcu- 
lation as in the previous section, as are the fluxes F^ +1 . The notation in Eq. 

(54) indicates that primitive variables are converted into conservative vari- 
ables which are then updated through conservative fluxes and then converted 
back into primitive form. Thus, if we indicate with AFj £ the undivided flux 
difference in the d direction for the total energy, the above transformations 
imply the following correction 



Ai^ -> AFy pE + {a-l)\ (p y+ . + PiJ ._i) (« s>w+ i 



This modification leads to a pressure change in accord to Eq. (41). Similarly, 
the entropy flux difference is modified as 



Af* - AF£ + (a - 1) ( 7 - 1) \ [{ P s\ j+ , + (ps).^ 



4 Stability Considerations 



The method outlined above satisfy a number of conditions required for nu- 
merical stability. It is easy to see from Eq. (40) that, as dA/de — > — oo, the 
internal energy decays rapidly to its equilibrium value, and thereafter remains 
constant, at that value. Inspection of the characteristic analysis shows that, 
in this limit, no information is carried along the entropy wave corresponding 
to the eigenvalue A . This means that the system of equations (1) effectively 
reduces to the equilibrium system in which the internal energy is fixed at its 
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equilibrium value. In addition, Eq. ( 36) and (37) indicate that the so called 
subcharacteristic condition for the characteristic speeds at equilibrium is al- 
ways satisfied. That is 

A_ < A c _ ff < A < Af < A+ (55) 

where and A +; o,- are the equilibrium and frozen eigenvalues, respectively. 
The above condition, while being necessary for the stability of our linearized 
system [19], also guarantees that the numerical solution tends to the solution 
of the equilibrium equation as the relaxation time tends to zero [3]. Since 
the structure of the equations and the numerical framework, including the 
Riemann solver, remain basically unaltered with respect to classic Godunov's 
schemes except for the modification of the sound speed, one expects the usual 
stability analysis to apply. The latter implies the familiar CFL condition on 
the time step 

max(|A*|)^ < 1 * = -0,+. (56) 
Ax 

As for the step involving the source update, stability analysis for deferred 
correction methods of the type adopted here was carried out through numerical 
experiments by Dutt et al. [5]. Minion [12] extends such considerations to the 
case of semi-implicit schemes as the one adopted here. While the stability 
and convergence properties of such schemes have not been fully elucidated 
analytically, the analysis of these authors suggest that they are in general very 
satisfactory and competitive with commonly employed modern integration 
schemes. 

Here we show that, provided that the CFL condition in Eq. (56) is satisfied, 
our method is A-stable, in the sense further specified below. To demonstrate 
this we apply the method to the following model problem [2] 

d ^l = AY + BY 
dt 

Y(0) = 1 

where Y : R — > C and A, B G C and represent the non-stiff and stiff part of 
the equation, respectively. Using the notation 

Y n+1 = P( Zl ,z 2 )Y n , 

where P(z 1 ,z 2 ) is the operator corresponding to the proposed method, Z\ = 
AAt, z 2 = BAt, the stability region of the method, P, is defined as the region 
Sp = {zi,z 2 G C : \P(zi,z 2 )\ < 1}. A method, P, is A-stable if Sp includes 
the plane C_ = {z G C : R(z) < 0}. Inspection of equations (9)-(12) and 
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simple algebraic manipulation lead to the expression 



P(z 1 ,z 2 ) = [l + P h (z 1 )] 



, 3 z 2 

1 Z2 — 

2 , 2 



1 - z 2 



(57) 



where is the hydro dynamic operator given by Godunov's method. Using 
the CFL conditions, which ensures |P(zi,0)| = |1 + Ph(zi)\ < 1, we find: 



\P(Z!,Z 2 )\ 2 < 



(i - z 2 y 



< 1 V^ 2 : R(z 2 ) < 0. 



(58) 



5 Extension to the Case A p ^ 



When the source term depends on both the internal energy and the gas density, 
A p 7^ and we obtain 



1 



a -I) 





1 

a 



(59) 



with a defined in Eq. (33). As a result 



A eS = 





I (9P) 
P \9pJ ( 







I (dp) 
P \9eJ l 

+ 







+ ul 



(60) 



and the sound speed is now given by 

A 



Ceff = 



P 



A e p 




(61) 



Since a < 1 the term in squared brackets can become negative and the sound 
speed imaginary. This behavior is related to the fact that when A p ^ the gas 
is prone to thermal instability so that the scheme cannot be simply generalized 
without taking into account the specific properties of the source term. In 
general one cannot expect an implicit method to work properly except in the 
case of a system with a stable solution. For a 7— law equation of state, Eq. (38), 
c 2 cS > requires 

e A e 1 — a , . 

PA^^T) (62) 
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which is reminiscent of the thermal stability criterion [6], in which case the 
term on the right-hand-side is 1. In both the stiff limit and non-stiff limits the 
RHS in Eq. (62) is of order — A e At >> 1, indicating the potential for triggering 
thermal instability of 'numerical nature'. For example, consider a source of the 
form A(p, e) = p n A(p, e), so that 

A p (p, e) = np- 1 A(p, e) + p n A p (p, e). (63) 

In general the former term can take both positive and negative values. So 
even though it vanishes at equilibrium, its effect is destabilizing and should be 
resolved in time. Depending on the definition of A, it is possible that A p > 0. 
Only in this case is the latter term stabilizing and should contribute to the 
sound speed in Eq. (61). 

So our approach is to decouple any destabilizing component of A p , which we 
indicate with A Pi< , from the characteristic analysis and associate it explicitly 
with the source term so that its effect does not enter the sound speed. In this 
case one would have to add a term 

Ap = pAe = -(a - 1) ^ p 2 (V • u) ^ (64) 

to the pressure component of the right hand side of Eq. (52). In order to 
preserve second order accuracy, one would require that the above term is 
resolved in time, i.e. the time step is sufficiently small that Ae < e. The 
restriction placed by Eq. (52) depends on the shape of the source function. 
However, near equilibrium it does not play any role because by definition to 
zeroth order the source term is zero. In fact we find that in all test cases with 
a density dependent source explored below, including the one in section 6.3, 
the condition (52) never constrained the time step. 



6 Tests 



In this section we test the performance of the proposed method in terms of 
both accuracy and robustness. As for the accuracy we consider a set of one 
dimensional problems for which the analytic solution is known. In particu- 
lar, we use the test problems in [14] for an isothermal rarefaction fan and an 
isothermal shock wave and consider a flow with a stiff relaxation term in the 
limit in which the relaxation time approaches zero. We then consider the case 
of a sinusoidal perturbation with wave-vector both parallel (1-D) and skew 
(2-D) with respect to the x-axis, and prove the second order accuracy of the 
scheme for a variety of stiffness conditions. This we show both for the case in 
which the source term does or does not depend on density. As for the robust- 
ness of the method, we turn to multidimensional problems involving strong 
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shocks and large spatial gradients. In particular we consider the interaction 
of a strong shock with a spherical cloud, again assuming a variety of stiffness 
conditions. The aim of the tests is to prove the code performance in the case 
of complex and computationally more challenging calculations. 

As for the source term, in the following we mostly present results based on a 
relaxation law of the form 



A = -K p c 




(65) 



where K is the heat transfer coefficient and the internal energy, e, is related to 
pressure and density by the equation of state (38), and ( and 77 are parameters. 
When ( = r] = the relaxation law expressed by Eq. (65) reduces to the case 
studied in [14] and in the limit K — > 00 it enforces isothermality. 

We test the case of a density dependent source term, by setting either ( or 
77, or both parameters, to a non zero value. Only the latter case is reported 
here although in all cases we obtain consistent results in terms of convergence 
and accuracy. When rj 7^ 0, Eq. (65) forces the system towards an equilib- 
rium configuration described by polytropic-like equation of state in which 
e = e (p/p ) v . Thus, when r\ > 0, Eq. (61) implies an effective adiabatic 
index that, as it should be, tends to (1 + rj), as a — > 0. 



6.1 Riemann Problems 



We first consider one-dimensional Riemann problems described by the follow- 
ing initial conditions 



(p,u,p)[x,t = 0] 



(pi,u h pi) ifx<0.5 
(p r , u r , p r ) if X > 0.5 



(66) 



and with a source term described by Eq. (65). Following [14] we adopt 



e 



K = W 

Pl,r 



Pl,r (7-1) 

7 = 1.4 
Ax = 2.5 x 10" 3 . 



(67) 
(68) 

(69) 
(70) 



The stiff nature of the problem is apparent as K 1 <^ Ax/c c s, i.e. the relax- 
ation time is much shorter than the hydrodynamic time scale. 
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Fig. 1. Isothermal rarefaction wave. From top to bottom: density, velocity and pres- 
sure solutions, respectively. Open dots and solid line indicate the numerical and 
analytic solution, respectively. The initial conditions are given in Eq. (71) with 
ui = —0.8. A mesh size Ax = 2.5 x 10~ 3 was employed. 

6.1.1 Isothermal Rarefaction 

We begin by setting the state variables to the values 
Pl = 1.0, Pl = 0.4, m = -0.8 

HI HI , (71) 

p r = 2.5, p r = 1.0, u T = ui + 0.5795 

representing an isothermal rarefaction in the A + characteristic family. For the 
calculation we employ a grid with N ce u = 400 grid cells [14]. The results from 
the code (open dots) are illustrated together with the analytic solution (solid 
line) in Fig. 1. From top to bottom the plot shows the density, velocity and 
pressure solutions at time t = 0.4, respectively (the same time as in [14]). 
The solutions are free of numerical artifact and well reproduce the analytic 
solution. In particular both the foot and the front edge of the rarefaction 
wave are accurately reproduced as sharp features. In addition, there is no 
numerical 'kink' along the wave in correspondence of the sonic point, that 
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is the the eigenvalues A + = u + c e s changes sign 1 as it was noticed in the 
'non-stiff' schemes presented for comparison in [14]. If we estimate the error 
in the numerical solution as in [14] 



N, 



cell 



N 



(72) 



that is the average of the absolute value of the difference between the numerical 
and analytic result, we find that the error is e p = 4.2 x 10~ 4 for the density 
and e u — 1.5 x 1CT 4 for the velocity. The latter is a factor almost 20 smaller 
than obtained with the 'frozen method' proposed in [14], most likely owing to 
the sharper resolution of our method at the rarefaction front and foot. This 
is visible from comparing the analytic and numerical solutions in Fig. 1. It 
is also consistent with the norm of the errors (see Eq. (79) in Sec. 6.2), 
which gives Halloo = 7.3 x 10~ 3 and ||£ p ||oo = 1-6 x 1CT 2 , indicating a localized 
error as opposed to one that is uniformly distributed. 



6.1.2 Shocks 

Next we study the case of an isothermal shock with initial conditions 

p T = 1.0, p r = 0.4, u r = u r 

pi = 2.5, pi = 1.0, ui = u r + 0.6. 

The numerical results (solid dots) are shown together with the analytic solu- 
tion (solid line) in Fig. 2 for two different values of u r , namely —1.2 (top left) 
—0.3 (top right) producing isothermal shock fronts slowly moving to the left 
and the right, respectively. Neither artificial viscosity nor flattening was em- 
ployed and a van Leer's limiter was used. Overall the algorithm performs very 
well. The shock positions accord with the analytic value. The shocks are well 
captured within a couple of zones, indicating that the properties of the scheme 
have not degraded with respect to the non-stiff case. We notice that minor 
oscillations appear in a few zones downstream the left moving shock front. 
These have not been introduced by our method for treating the stiff source, 
but rather, are due to the fact that the dissipation in a Godunov method 
vanishes for slowly-left-moving shocks, such as the one being computed here. 
A thorough discussion on this is found in [20]. In particular, we find that the 
same oscillations appear in the purely hydrodynamic version of the algorithm 
with the relaxation term turned off, if we use an adiabatic index < 7 — 1 1 
in order to mimic isothermality. 



This occurs as the effective sound speed is c c g ~ 0.63 and ui varies from —0.8 to 
-0.2205. 
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Fig. 2. Top panels: slow left-moving (left panel) and right-moving (right panel) 
isothermal shock waves. In each panel, from top to bottom: density, velocity and 
pressure solutions, respectively. Filled dots and solid line indicate the numerical 
and analytic solution, respectively. The initial conditions are given in Eq. (73) 
with u r = —1.2 (top left), and u r = —0.3 (top right). Bottom: slow right-moving 
quasi-isothermal shocks. In both cases we set u r = —0.3 and use a heat transfer 
coefficient K = 500 (bottom left) and K = 50 (bottom right) respectively. The 
solid line in the bottom right panel corresponds to a solution obtained with a much 
higher resolution run using Ax = 8.3 x 10~ 5 so that the reaction time is resolved. 
In all other cases a mesh size Ax = 2.5 x 10~ 3 was employed. 



Finally, in the bottom right panel the same initial conditions as in the top 
right panel are used but in combination with a much smaller heat transfer 
coefficient, K = 500 (bottom left) and K = 50 (bottom right). In these 
cases K^ 1 < Ax/c e g and K^ 1 > Ax/c c s, respectively, so that while the 
gas behavior is not strictly isothermal, the relaxation time is still relatively 
short. Assessing the algorithm performance for this situation is of relevance 
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as well, as in general stiffness of the conditions will vary across a flow. The 
bottom panels of Fig. 2, in addition to density, velocity and pressure also 
present results for the temperature. As it appears from this plot, the numerical 
solution is very satisfactory, without numerical artifacts or oscillations. For 
the case K = 50, we take the test one step further and compare the result 
obtained using the current grid settings (solid dots), that is N ce u = 400, Ax = 
2.5 x 10~ 3 , with those from a much higher resolution run (solid line) in which 
N ce n = 12000, Ax = 8.33 x 10~ 5 so that the reaction time is fully resolved. 
This comparison nicely shows that our modified scheme is able to capture 
the correct behavior of the flow even at intermediate stiffness conditions away 
from a purely isothermal behavior. 

6.2 Convergence Rates in Smooth Flows 

In this section we test the convergence of the method by studying the case of 
a smooth flow with the following initial conditions 



where r is the position vector and we use p = 7 = 1.4. The above initial 
conditions produce a sinusoidal wave with amplitude A propagating in the 
domain along the direction defined by the vector k. While we have experi- 
mented with various values for the parameters A, k and K, below we present 
results for a few cases only, summarized in Table 1. 

In particular we consider a perturbation amplitude A = 10~ 2 and both a 
wave-vector aligned with the grid k = (1,0) and skew with respect to it, 
k = (2/y/E, l/y/h). We adopt a source term as given in Eq. (65) with values 
of the transfer coefficient K — 1, 50, 10 8 to explore different regimes in which 
the relaxation is resolved in time, is stiff as well as the intermediate regime 
(cases A-F). We then repeat case C and F but with the initial value of the 
internal energy offset from the equilibrium value by Se/eo = 40% (cases G-H). 
Finally, we consider the case in which the source term depends both on density 
and internal energy, as described by Eq. (65). In particular, we show results 
concerning the case in which K = 10 8 , £ = 1, rj — 0.1, 5e/e = 40% (cases 
I-L). Consistent test results were also found by setting ( — 1, rj — as well as 



In order to measure the rate at which the numerical solution converges, for 
each problem we carry out a set of 5 simulation runs employing N ce u = 



P = Po + — [cos(2tt k • r) + 1] 

P = Po = 0.5 
Ux = UxO = 0.3 

Uy = UyO = 0.5 



(75) 
(76) 
(77) 



(74) 



C = 0,?7 = 0.1. 
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Table 1 

Run Set for Convergence Study with Relaxation Law Eq. (65) 



run 


A 


k 


if 


C 


V 


Note 


A 


1(T 2 


(1,0) 


if = 1 










B 


lO- 2 


(1,0) 


if = 50 










C 


KT 2 


(1,0) 


if = 10 8 










D 


io- 2 


{2/y/E, 


if = 1 










E 


KT 2 


(2/V5, 1/a/5) 


if = 50 










F 


KT 2 


(2/V5, 1/V5) 


if = 10 8 










G 


io- 2 


(1,0) 


if = 10 8 








5e/e = 0.4 


H 


1(T 2 


(2/V5, 1/V5) 


if = 10 8 








5e/e = 0.4 


I 


io- 2 


(1,0) 


if = 10 8 


1 


0.1 


5e/e = 0.4 


L 


io- 2 


(2/V5, 1/V5) 


if = 10 8 


1 


0.1 


5e/e = 0.4 



32, 64, 128, 256, 512 for a total range of 32. Note that the stiffness conditions 
do not change significantly as the grid is refined within the range of resolutions 
considered here. Also, the smallness of the perturbations is such that the term 
given by Eq. (64), to be added to the energy in the predictor step, is resolved 
in time 

The convergence rate is measured using Richardson extrapolation. Given the 
numerical result q r at a given resolution r we first estimate the error at a given 
point as 

£ r ;i,j = q r (hj) - q r +i(hj) ( 78 ) 
where q r +\ is the solution at the next finer resolution, properly spatially aver- 
aged onto the coarser grid. We then take the n-norm of the error 

L-n H^rlln ( ^ ] \^r;i ,j I ) C^9) 

where, v it j = Ax 2 is the cell volume, and estimate the convergence rate as 

Rn ~ \n(Ax r /Ax s ) ■ (80) 

For each studied case listed Table 1, we produce a corresponding table (2-5) 
reporting the Li, L 2 and norms of the error as defined above. Inspection of 
their values shows that the error drops with second order accuracy, supporting 
our analysis in Sec. 2.. 

A final experiment is designed to further prove that in the stiff limit (if = 
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Table 2 

Convergence Rates: 1-D Case: A = 10~ 2 , k = (1,0) 



density momentum 



iVcell 


Li 


L 2 


Loo 


R\ 


Li 


L 2 


Loo 


R\ 










K = 


1 








32 


4.3E-7 


9.5E-7 


2.8E-6 




6.3E-7 


1.4E-7 


4.0E-6 




64 


1.1E-7 


2.4E-7 


7.0E-7 


2.0 


1.3E-7 


2.8E-7 


8.0E-7 


2.3 


128 


2.7E-8 


6.0E-8 


1.8E-7 


2.0 


2.8E-8 


6.3E-8 


1.8E-7 


2.2 


256 


6.8E-9 


1.5E-8 


4.4E-8 


2.0 


6.7E-9 


1.5E-8 


4.2E-8 


2.1 



K = 50 



32 


4.0E-7 


8.8E-7 


2.6E-6 




7.8E-7 


1.7E-6 


4.9E-6 




64 


1.1E-7 


2.5E-7 


7.2E-7 


1.9 


1.5E-7 


3.2E-7 


9.2E-7 


2.4 


128 


3.1E-8 


6.9E-8 


2.0E-7 


1.8 


2.9E-8 


6.5E-8 


1.8E-8 


2.4 


256 


8.5E-9 


1.9E-8 


5.4E-8 


1.9 


6.2E-9 


1.4E-8 


3.9E-8 


2.2 










K = 


10 8 








32 


4.1E-7 


9.2E-7 


2.7E-6 




8.2E-7 


1.8E-6 


5.91-6 




64 


9.7E-8 


2.1E-7 


6.4E-7 


2.1 


1.6E-7 


3.6E-7 


1.0E-6 


2.4 


128 


2.3E-8 


5.2E-8 


1.5E-7 


2.1 


3.6E-8 


8.0E-8 


2.3E-7 


2.1 


256 


5.7E-9 


1.3E-8 


3.8E-8 


2.0 


8.5E-9 


1.9E-8 


5.4E-8 


2.1 



i?i is the convergence rate based on the L\ errors. 



10 8 ) the proposed scheme converges to the correct asymptotic (isothermal) 
behavior. To do that we employ a Godunov method for the isothermal fluid 
equations and run again case C in Table 1. A comparison of the solutions 
obtained with the isothermal code and our proposed method is reported in 
Table 6. It shows that the difference between the two is always negligible with 
respect to the estimated truncation error, thus validating our convergence 
study. 
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Table 3 

Convergence Rates: 2-D Case: A = 10~ 2 , k = {2/y/E, 1/yg) 



density momentum 



iVcell 


Li 


L 2 


Loo 


R\ 


Li 


L 2 


Loo 


R\ 










K = 


1 








32 


3.2E-5 


3.5E-5 


5.3E-5 




4.3E-5 


4.8E-5 


7.0E-5 




64 


8.6E-6 


9.5E-6 


1.4E-5 


1.9 


9.8E-6 


1.1E-5 


1.6E-5 


2.1 


128 


2.2E-6 


2.5E-6 


3.6E-6 


2.0 


2.3E-6 


2.6E-6 


3.8E-6 


2.1 


256 


5.7E-7 


6.4E-7 


9.2E-7 


2.0 


5.7E-7 


6.3E-7 


9.1E-7 


2.0 



K = 50 



32 


6.2E-5 


7.0E-5 


1.0E-5 




3.6E-5 


4.0E-5 


5.9E-5 




64 


1.2E-5 


1.4E-5 


2.0E-5 


2.4 


7.9E-6 


8.8E-6 


1.3E-5 


2.4 


128 


2.7E-6 


3.0E-6 


4.3E-6 


2.2 


1.8E-6 


1.9E-6 


2.8E-6 


2.1 


256 


6.2E-7 


6.9E-7 


1.0E-6 


2.1 


4.0E-7 


4.5E-7 


6.5E-6 


2.2 










K = 


= 10 8 








32 


7.4E-5 


8.2E-5 


1.2E-4 




4.5E-5 


5.0E-5 


7.3E-5 




64 


1.6E-5 


1.8E-5 


2.6E-5 


2.2 


1.0E-5 


1.1E-5 


1.7E-5 


2.2 


128 


3.8E-6 


4.2E-6 


6.1E-6 


2.1 


2.5E-6 


2.8E-6 


4.1E-6 


2.0 


256 


9.4E-7 


1.0E-6 


1.5E-6 


2.0 


6.2E-7 


6.9E-7 


9.9E-6 


2.0 


Ri 


is the convergence rate based on 


the L\ 


errors. 









6.3 Adaptive Mesh Refinement and Strong Shock Problems 

In applications involving the interaction of strong shocks, it is useful to use the 
dissipation mechanisms described in [4], which generalize without modification 
to the present case. In addition, it is also desirable to couple this method to 
a block-structured adaptive mesh refinement (AMR) [1,10]. In AMR calcula- 
tions, the conservative variables are updated for the conservative fluxes in two 
steps. The first step constitutes the main flux update and it simply consists 
in modifying the state variables U for the total fluxes across the cell bound- 
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Table 4 

Convergence Rates: Off Equilibrium Case: 8e/eo = 0.4, A = 10~ 2 , K = 10 8 



density momentum 



iVcell 


Li 


L 2 


Loo 


R\ 


Li 


L 2 


Loo 


R\ 










k = 


(1,0) 








32 


4.8E-7 


1.1E-6 


3.1E-6 




7.6E-6 


1.7E-6 


4.8E-6 




64 


1.0E-7 


2.3E-7 


6.9E-7 


2.2 


1.5E-7 


3.3E-7 


9.5E-6 


2.3 


128 


2.4E-8 


5.4E-8 


1.6E-7 


2.1 


3.3E-8 


7.3E-7 


2.1E-7 


2.2 


256 


5.9E-8 


1.3E-8 


3.9E-7 


2.0 


7.7E-9 


1.7E-8 


5.9E-8 


2.1 


k = (2/v% 1/V5) 


32 


7.4E-5 


8.2E-5 


1.2E-4 




4.4E-5 


4.9E-5 


7.2E-5 




64 


1.6E-5 


1.8E-6 


1.7E-5 


2.2 


1.0E-5 


1.1E-5 


1.7E-5 


2.1 


128 


3.8E-6 


4.2E-6 


6.2E-6 


2.1 


2.5E-6 


2.8E-6 


4.0E-6 


2.0 


256 


9.4E-7 


1.0E-6 


1.5E-6 


2.0 


6.1E-7 


6.8E-7 


9.9E-7 


2.0 



t i?i is the convergence rate based on the L\ errors. 



aries. In addition, as part of the operations of synchronization among different 
levels, the conservative variables at the coarse-fine grid interfaces are further 
updated for the flux difference between the level on which they are defined 
and the next finer level. This operation is referred to as refluxing and it is 
aimed at preserving the conservative character of the numerical scheme when 
applied to a hierarchy of nested grids. 

For the purpose of the current discussion, the effect of this operation can be 
expressed as 

U - ^5F (81) 
Ax v ' 

where 5F is the difference between the fluxes at the coarse-fine interface com- 
puted on a given level and the next finer level. In AMR calculations refluxing 
on a given level is enforced as a separate operation, after the source update 
and the main flux update have been carried out on that level and also on all 
finer levels. Therefore, an additional measure must be taken to ensure that the 
effects of refluxing are also subjected to the action of the deferred corrections 
(just like the flux update does). Thus, inspection of Eq. (9)-(ll) indicates that 
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Table 5 

Convergence Rates: Off Equilibrium, p-dependent Source Case: 5e/eo = 0.4, ( = 
l,i] = 0.1, A = 10~ 2 , K = 10 s 



density momentum 



Wcell 


Li 


L 2 


Loo 


R\ Li 


L 2 


Loo 




k=(l,0) 


32 


2.6E-7 


5.7E-7 


1.8E-6 


9.0E-7 


2.0E-6 


5.7E-6 




64 


6.1E-8 


1.4E-7 


4.2E-7 


2.1 1.9E-7 


4.1E-7 


1.2E-6 


2.2 


128 


1.5E-8 


3.3E-8 


1.0E-7 


2.0 4.2E-8 


9.3E-8 


2.7E-7 


2.2 


256 


3.5E-9 


7.7E-9 


2.4E-8 


2.1 l.OE-8 


2.2E-8 


6.5E-8 


2.1 


k = (2/^5, 1/V5) 


32 


6.6E-5 


7.4E-5 


1.1E-4 


5.1E-5 


5.7E-5 


8.3E-5 




64 


1.5E-5 


1.7E-5 


2.5E-5 


2.1 1.2E-5 


1.3E-5 


1.9E-5 


2.1 


128 


3.7E-6 


4.1E-6 


5.9E-6 


2.0 2.8E-6 


3.1E-6 


4.6E-6 


2.1 


256 


9.1E-7 


1.0E-6 


1.5E-6 


2.0 7.0E-7 


7.7E-7 


1.1E-6 


2.0 


R\ is 


the convergence rate based on the L\ errors. 








Table 6 

Comparison with 


a Purely Isothermal Solution: A = 10 2 


, k=(l,0) 








density 




momentum 




iVceU 


Li 


L 2 




Loo L\ 


L 2 




Loo 



32 


2.7E-11 


6.0E-11 


1.7E-10 


1.8E-10 


4.0E-11 


1.1E-10 


128 


2.7E-11 


6.1E-11 


1.7E-10 


1.8E-10 


4.0E-11 


1.1E-10 


256 


2.7E-11 


6.1E-11 


1.7E-10 


1.8E-10 


4.0E-11 


1.1E-10 


512 


2.7E-11 


6.1E-11 


1.7E-10 


1.8E-10 


4.0E-11 


1.1E-10 
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Fig. 3. Logarithmic pressure maps from the shock-cloud interaction run. The shock 
Mach number is 10 and the cloud overdensity is 10. The left panel shows the 
'isothermal' case with K = 10 8 and the right panel shows the case in which 
K~ l ~ Ax/v s h oc k, i.e. the relaxation time is comparable to the shock cell cross- 
ing time. These calculations were performed with an AMR code which employed a 
base grid of 256 x 256 zones and two additional levels of refinement with refinement 
ratio 2. 

the flux correction must be modified according to 

SF -> {(I - AtVuSlu,)- 1 + (I - AtVuSlu)- 1 [i - (I - AtV^S^)" 1 ] } 5F. 

(82) 

In the following we employ an AMR code and carry out a calculation involving 
the interaction of a spherical overdense region with a strong hydrodynamic 
shock to assess the robustness of our proposed numerical method. We assume 
a cloud overdensity with respect to the ambient medium x = 10 and a shock 
Mach number Ai = 10. We use a base grid of 256x256 zones and allow for 
two additional levels of refinement with refinement ratio 2 in regions where 
the undivided, relative density gradients, Ap/p, exceed 20%. 

We begin assuming that a stiff relaxation term of the form in Eq. (65) acts on 
the flow internal energy and we consider both the case of a exceedingly large 
transfer coefficient, K = 10 8 , as well as the case in which the relaxation time 
is comparable to the shock cell crossing time. This requires, roughly, that 

K~ x ~ Ax/u shock (83) 

where Ax is the mesh size and u s h oc k is the shock speed. Note that the stiffness 
is sufficiently large that refining by a factor of 2 or 4 does not make the 
problem significantly less stiff. At simulation start the temperature is constant 
throughout the domain, so that the cloud is in thermal equilibrium but it has 
higher pressure than its surroundings. As a result, it expands sonically into the 
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Fig. 4. Logarithmic density (left) and pressure (right) maps from the shock-cloud 
interaction run for a density dependent source term at t = 0.018 (top) and t = 0.07 
(bottom) time units. As before, the shock Mach number is 10, the cloud overdensity 
is 10 and the calculation was performed with an AMR code employing a base grid 
of 256x256 zones and two additional levels of refinement with refinement ratio 2. 

background. The shock propagates from the right to the left along the x-axis 
and as it runs into the cloud it crushes it. In Fig. 3 we plot the logarithmic 
pressure map as the shock is roughly half-way through the cloud. The high 
pressure postshock region is clearly thinner in the case of the larger value of 
K and the shock has also propagated slightly further down the axis. In both 
cases, however, the result is sound and shows no sign of numerical artifact 
both in the presence of strong shock and large gradients, and independently 
of the magnitude of the transfer coefficient. 

As a final test, we consider the same shock cloud interaction problem as de- 
scribed above but with a source function appropriate for a mixture of hydrogen 
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(76 %) and helium (24 %) illuminated by a uniform ionizing background ra- 
diation field. The cooling part of the source function is proportional to the 
density and the equilibrium temperature, of order 10 4 K, depends slightly on 
the density. This function has a very strong temperature gradient about the 
equilibrium value, behaving analogously to the stiff source terms used in the 
previous sections were accuracy and convergence studies were carried out. 

We set the background gas temperature to 10 6 K and its number density to 
0.1 cm~ 3 . The gas is collisionally ionized, its sound speed is of order 1.2 x 10 7 
cm s~ x and it has a cooling time r coo i = P/(7 — l)pA ~ 1.4 x 10 7 yr. With 
a box size L = 1.7 kpc = 5.2 x 10 21 cm, the latter is much longer than the 
CFL time, tcfl < Ax/u s h oc k — 5.4 x 10 3 yr. The cloud of overdense gas is in 
pressure equilibrium with a density contrast x = 10 so that its temperature is 
10 5 K. When unperturbed, the cloud's gas cooling time is about 1.4 x 10 4 yr 
~ 2.6 x t C fl- A background radiation field, producing about T ~ 2.4 x 10~ 12 
s -1 ionizations of neutral atoms of hydrogen and helium keeps the cloud's 
temperature at the equilibrium value of ~ 1.510 4 K. However, the cloud's 
pressure quickly falls below the background value and, as a minor effect, the 
cloud slowly contracts. 

Fig. 4, shows a snapshot of the density (left) and the pressure (right) dur- 
ing the initial (top) and final stages (bottom) of the simulation. The reverse 
shock is non-radiative, thus extending further ahead of the cloud than in the 
previous cases in Fig. 3. Inside the cloud strong radiative losses prevent the 
full temperature rise in the postshock region and produce a density jump sub- 
stantially larger than in the corresponding adiabatic case. The bottom panel 
shows the later stages of the cloud evolution, when Rayleigh- Taylor instability 
with scales comparable to the cloud size have developed and are shredding the 
cloud. As in the previous case, in which the source term is described by a re- 
laxation law, the code appears to produce reliable numerical results, without 
numerical artifact despite the presence of strong shock and large gradients. 



7 Conclusions 

We have presented a second order accurate semi-implicit predictor-corrector 
scheme to treat stiff source terms within the framework of higher order Go- 
dunov's methods. Our treatment of the predictor step for computing the hy- 
perbolic fluxes, is based on the derivation of a local effective dynamics us- 
ing Duhamel's formula. This leads to a conventional second-order Godunov 
method when the system relaxation time is larger than the time step and to 
a second-order Godunov method for the isothermal equations in the limit of a 
stiff source term. Finally, we obtain a semi implicit corrector using a one-step 
second-order accurate deferred corrections method as suggested in [5,12]. 
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Our tests indicate that the proposed method is stable, robust and its second 
order accuracy preserved across a variety of stiffness conditions. We have also 
discussed the case of a general source term which depends both on e and p 
and shown that the method is applicable provided that the flow is thermally 
stable or the non-stiff part of the source term is resolved in time. 

The additional cost involved in the formulation of our scheme is minimal; all 
it requires is an estimate of the term A e which in a purely relaxation case is 
trivial and for a more complicated source term (such as the case of radiative 
losses) is still minor compared to the estimate of the source term itself. In our 
implementation the factor (2(7 — 1) + 1 is stored as an additional primitive 
variable and used as polytropic index in the characteristic analysis in stead of 

7- 
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